Wang-Landau sampling for quantum systems: 
algorithms to overcome tunneling problems and calculate the free energy 
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We present a generalization of the classical Wang-Landau algorithm [Phys. Rev. Lett. 86, 2050 
(2001)] to quantum systems. The algorithm proceeds by stochastically evaluating the coefficients 
of a high temperature series expansion or a finite temperature perturbation expansion to arbitrary 
order. Similar to their classical counterpart, the algorithms are efficient at thermal and quantum 
phase transitions, greatly reducing the tunneling problem at first order phase transitions, and allow 
the direct calculation of the free energy and entropy. 
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Monte Carlo simulations in statistical physics now have 
a history of nearly half a century starting with the sem- 
inal work of Metropolis While the Metropolis algo- 
rithm has been established as the standard algorithm for 
importance sampling it suffers from two problems: the 
inability to directly calculate the partition function, free 
energy or entropy, and critical slowing down near phase 
transitions and in disordered systems. 

In a standard Monte Carlo algorithm a series of config- 
urations is generated according to a given distribution, 
usually the Boltzmann distribution in classical simula- 
tions. While this allows the calculation of thermal aver- 
ages, it does not give the partition function, nor the free 
energy. They can only be obtained with limited accuracy 
as a temperature integral of the specific heat, or by using 
maximum entropy methods . 

The problem of critical slowing down has been over- 
come for second order phase transitions by cluster up- 
date schemes for classical ^ and quantum systems 
[@7 l^l 0- For first order phase transitions and sys- 
tems with rough free energy landscapes a decisive im- 
provement was achieved recently by a new algorithm for 
classical systems due to Wang and Landau In con- 
trast to related methods - such as the multicanonical ||] 
or the broad histogram |^ method - this new algorithm 
scales to large systems, does not suffer from systematic 
errors and needs no a priori knowledge. The key idea 
is to calculate the density of states p{E) directly by a 
random walk in energy space instead of performing a 
canonical simulation at a fixed temperature. By visit- 
ing each energy level E with a probability l/p{E), this 
algorithm achieves a fiat histogram and good precision 
over the whole energy range. Besides being efficient at 
first and second order phase transitions this algorithm 
allows the direct calculation of the free energy from the 
partition function Z = p{E)e~^^''^'^ . The internal 
energy, entropy, specific heat and other thermal proper- 
ties are easily obtained as well, by diflFerentiating the free 
energy. Within a year of publication this algorithm has 



been improved using A^-fold way and multibondic 
|l2| sampling schemes, has been applied to Potts mod- 
els [|l^ and generalized to reaction coordinates jlj] , con- 
tinuum models |]l5| , polymer films |l^ ], and to protein 
folding 0. 

Since simulations of quantum systems suffer from the 
same problems as classical simulations, in particular from 
long tunneling times at first order phase transitions and 
the inability to calculate the free energy directly, an ex- 
tension of this algorithm to quantum systems is highly 
desired. Here we present two such algorithms. The first 
one is based upon a high temperature series expansion. 
Similarly to the classical algorithm it allows the calcu- 
lation of the free energy as a function of temperature, 
making it ideal for the study of thermal phase transitions. 
The second algorithm renders the high temperature se- 
ries expansion into a perturbation expansion in one of the 
coupling constants and is suitable for the investigation of 
quantum phase transitions. 

Quantum Monte Carlo algorithms start by mapping 
the quantum system to a classical system. This can be 
done either through a discrete or continuous time j|] 
path integral representation or by a stochastic series ex- 
pansion (SSE) in the inverse temperature [|l9). While our 
algorithms can be formulated in both representations, we 
here present the SSE version as it is the easiest and most 
natural representation for most problems. 

We start by expressing the partition function as a high 
temperature expansion 

oo „„ CXD 

Z = Tve-f'" = J2 ^Tr(-i/)" ^ g(n)/3", (1) 

n=0 ' n=0 

where the n-th order series coefficient g{n) = 
Tr(— TJ)"/??,! will play the role of the density of states in 
the classical algorithm. The algorithm performs a ran- 
dom walk in the space of series expansion coefficients, 
achieves a fiat histogram in their orders n and calculates 
the coefficients g{n). Employing Eq. (|^) we can then cal- 
culate the free energy, internal energy, entropy and spe- 
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cific heat directly. Thermal averages of observables can 
be measured as in conventional Monte Carlo algorithms 
by recording a separate histogram for the expectation 
values at each order. 

Next we note that in a computer simulation the series 
expansion (|l]) needs to be truncated at an order A. Since 
the orders relevant for a given inverse temperature (3 are 
sharply peaked around P\E{P)\, where E{(3) is the mean 
energy at inverse temperature /3, this cutoff does not in- 
troduce a systematic error. Its main consequence is to 
restrict the accessible temperature range to /3 < A/E{(3). 

The next step is to introduce a complete set of ba- 
sis states {|a)}, and to decompose the Hamiltonian H 
into diagonal and offdiagonal bond terms Hj^'^"' . For sim- 
plicity we restrict the following discussion to a spin-1/2 
Heisenberg model where this decomposition for a bond 
b = reads i/g) = JS^S^ - J/4, and i/gj) = 

J/2{S+S- + SrS+). The offset -J/4 is added to the 
diagonal part in order to render the matrix elements of 
—H nonnegative. Using this decomposition, we can write 
the partition function as |l^ 

a {Sa} ' »=0 

where the operator string 5a = ((6i, ai), . . . , (^a, oa)) 
is a concatenation of n bond operators and A — n unit 
operators. 

Comparing Eq. (|l|) to Eq. (H) we see that we can ob- 
tain g{n) by counting the number of times a configuration 
with n non-unit operators is observed during a simula- 
tion at an inverse temperature /3 — 1. As the dynamic 
range of g(n) spans thousands of orders of magnitude 
[g{0)/g{A) is up to 10^'^°''° for the examples given below] 
simply collecting a histogram will not be effective. There- 
fore a variant of the classical Wang-Landau method 
will be employed: by reweighting a configuration of n-th 
order with l/g{n) a flat histogram of the orders n can be 
achieved, thus sampling all orders equally well. 

Since g{n) is initially unknown we start with the (bad) 
guess g(n) — 1. Each time a configuration of n-th order 
is visited, g{n) is multiphed by a factor /, i.e. g{n) ^ 
fg{n). In our implementation we store the logarithms of 
these quantities to avoid overflow problems. The random 
walk is performed until the histogram H{n) - counting 
the number of times a configuration with n operators 
is observed - is reasonably flat. Similar to the classical 
case a maximum deviation of 20% from the mean value 
turned out to be reasonable. The multiplicative increase 
of g(ri) is essential for the success of the algorithm. Only 
that way the large range of g{n) can be mapped out in 
reasonable time, and g{n) converges rapidly to a rough 
estimate of the true distribution. Once the histogram 
is flat, it is reset to zero, / is decreased, in our case 
by / <— -v/Ji ^iid the process starts again, refining g{n) 
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FIG. 1: Free energy F , entropy S and specific heat C of 
an A'' = 10 site antiferromagnetic Heisenberg chain. Solid 
lines correspond to the MC results, indistinguishable from the 
dotted lines for the exact results. Also shown is the relative 
error e{F) of F compared to the exact result. 



further with smaller steps. This procedure is repeated 
until / gets as small as exp(10~®), so that an accurate 
estimate of g{n) with only negligible systematic errors 
will be available. The accuracy of the free energy and 
other calculated quantities is given by the statistical error 
which, as usual, scales with I/VTVmc where A^mc is the 
number of Monte Carlo steps. The overall normalization 
of g{n) follows from g(0) being equal to the dimension of 
the Hilbert space, which for a spin-1/2 Heisenberg model 
on a lattice with N sites is 2^. The initial choice of / is 
very important. Too small starting values result in long 
computation times, while too large values give extreme 
fluctuations in the initial iterations. As in the classical 
Wang-Landau method a good choice is to let /^"c i-,g 
of the same order of magnitude as the total number of 
configurations, which is of order 2^N^. Since usually 
A > A^, a good initial value is In / « ( A In A^)/AfMC • 

To finish the description of the algorithm we discuss 
the update steps in more detail. Any of the known update 
algorithms, employing local jl^, or cluster Q updates 
can be used. The only change in the acceptance proba- 
bilities from standard SSE algorithms is to set /3 — 1 and 
to include an additional factor g{n)/g{n') in the accep- 
tance probability for any move that changes the number 
of operators from n to n'. As an example we discuss 
the Heisenberg antiferromagnet. There the optimal algo- 
rithm is the loop algorithm, which in the SSE represen- 
tation consists of two parts: diagonal updates and loop 
updates. In a diagonal update step the operator string 
positions I = 1,...,A are traversed in ascending order. 
Empty and diagonal operators can be exchanged with 
each another. Using the notation \a{l)) — Y[\=i ^b" '*!'^) 
for the state obtained by acting on |a) with the first I 
bond operators, and M for the total number of interact- 
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FIG. 2: Scaling plot of the staggered structure factor of a 
cubic antiferromagnet as a function of temperature, obtained 
from simulations at a fixed temperature for various lattice 
sizes. The inset shows the specific heat as a function of tem- 
perature. The cutoff A = 500(L/4)^ restricts the accessible 
temperature range to T > 0.4J. 



ing bonds on the lattice, the update probabihties are 



g{n)M{c 



{d)\ 



g{n+l){k-n) 
g{n){A - n+1) 

g{n~l)M{a{l)\Hl'^\a{l)) 



(3) 



where P > 1 is interpreted as P = 1. This choice of 
update probabihties requires the least changes to an ex- 
isting SSE program. Alternatively the factors M, A — n 
and A — n + 1 can be dropped from the update probabil- 
ities, thus simplifying the algorithm. To correct for this 
omission, the obtained g{n) must then be multiplied by 
M"(A — n)!/A!. At each level I, independent of whether 
an update was performed or not (e.g. when the operator 
is off-diagonal) g{n) for the current value of n is incre- 
mented by /. The second part of the update cycle, the 
loop update, changes diagonal to off-diagonal bond op- 
erators without changing n and can be performed as in 
standard SSE algorithms. We refer to Ref. ^ for details. 

As a first example we show in Fig. |] results of calcu- 
lations for the free energy P, entropy S and specific heat 
C of an = 10 site antiferromagnetic Heisenberg chain, 
and compare to exact results. Using 10® sweeps, which 
can be performed in less than five hours on an 800 MHz 
Pentium-Ill CPU, the errors can be reduced down to the 
order of 10"''. The cutoff was set to A = 250, restricting 
the accessible temperatures to T > 0.05J. The sudden 
departure of the Monte Carlo data from the exact values 
below this temperature clearly shows this limit, which 
can be pushed lower by increasing A. The sudden devi- 
ation becomes even more pronounced in larger systems 
and provides a reliable indication for the range of validity 
of the results. 



FIG. 3: Average tunneling times (in units of Monte Carlo 
sweeps) between horizontal and vertical arrangement of 
stripes in a hard core boson model at a ratio V2/t = 3 of 
next nearest neighbor repulsion to kinetic energy on a 8 x 8 
sites lattice. The solid line is obtained using the standard SSE 
algorithm with directed loop updates. The dashed line is ob- 
tained using our algorithm, where the temperature is defined 
as the lowest temperature accessible in the simulation. 



To illustrate the efhciency of the algorithm close to 
a thermal second order phase transition, we consider 
the Heisenberg antiferromagnet on a simple cubic lat- 
tice. From simulations of systems with sites, L — 
4, 6, 8, 12, 16, we can calculate the staggered structure 
factor S'(7r,7r,7r) for any value of the temperature using 
the measured histograms. Figure ^ shows the scaling plot 
of S'(7r,7r,7r)/p2-'? with rj = 0.034. The estimate for the 
critical temperature Tc ~ 0.947 J, obtained in only a cou- 
ple of days on an 800 MHz Pentium-Ill CPU, compares 
well with earlier estimates pO{ . 

Next we demonstrate the efficiency of the algorithm 
at a first order phase transition by considering two- 
dimensional hardcore bosons with next-nearest neigh- 
bour interactions ||2^ . At low temperature and half filling 
this model is in an insulating phase with striped charge 
order and provides the simplest quantum mechanical 
model with stripes. We are currently investigating the 
thermal and quantum melting transitions of this stripe 
phase and probe for the existence of a nematic phase 
p2| . Simulations with conventional update schemes suf- 
fer from exponentially increasing tunneling times needed 
to change the stripe orientation from a horizontal to a 
vertical arrangement. The flat histogram in the order n 
in our algorithm reduces the tunneling times by many 
orders of magnitude already on small lattices (c.f. Fig. 
^) which demonstrates the efficiency of our algorithm at 
first order phase transitions. 

We now turn to our second algorithm, which applies to 
quantum phase transitions. Instead of scanning a tem- 
perature range we vary one of the interactions at fixed 
temperature. Defining the Hamiltonian as H = Hq + XV 
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FIG. 4: Scaling plot of the staggered structure factor of a 
Heisenberg bilayer as a function of the coupling ratio A = 
J/J'. Results are shown for various linear system sizes L. 
The temperature was chosen jSJ' = 2L, low enough to be 
in the scaling regime. The cutoff A = was chosen large 
enough to cover the coupling range J / J' < 1. The dynamical 
critical exponent of this model is z = 1 and r) = 0.034. 



well as thermodynamic averages accurately over a whole 
temperature or coupUng range. The algorithms can be 
used with any of the update schemes developed for the 
SSE representation and require only minor modifications 
to existing programs. Parallelization of the algorithms 
can be done like in the classical case by splitting the n- 
range into multiple random walks over a shorter range. 

Our algorithms are efficient not only at second or- 
der phase transitions but also at first order ones, where 
standard local and cluster update methods fail. The al- 
gorithms open up new possibilities for quantum Monte 
Carlo simulations: similar to the classical case we expect 
the algorithms to be efficient also for disordered systems 
and work is in progress to apply the methods to quantum 
spin glasses. Also, the ability to calculate the free energy 
and entropy will be useful for investigations of entropy- 
driven phase transitions, such as the reentrant melting 
transition observed in bosonic systems and anisotropic 
quantum magnets ||2^ . 

We thank D.P. Landau for discussions and acknowl- 
edge support of the Swiss National Science Foundation. 



we rewrite the partition function Eq. (Q) as 



an 



where on the right hand side we have collected all terms 
associated with A"* into g{nx). A similar algorithm can 
now be devised for this perturbation expansion up to 
arbitrary orders by setting A = 1, replacing M by j3M 
and S'(n) by g{n\) in Eq. (^. To normalize g{n\) there 
are two options. If i/o can be solved exactly, g{Q) can be 
determined directly. Otherwise, the normalization can 
be fixed using the first algorithm to calculate Z{(3) at 
any fixed value of A. Finally, even without normalization 
we can still obtain entropy and energy differences. 

We consider as an example the quantum phase tran- 
sition in a bilayer Heisenberg antiferromagnet whose 
ground state changes from quantum disordered to Neel 
ordered as the ratio A = J / J' of intra-plane (J) to 
inter-plane (J') coupling is increased ||2^. From the his- 
tograms generated within one simulation we can calcu- 
late the staggered structure factor 5(7r, tt) of the system 
at any value of A. In Fig. ^ we show a scaling plot of 
S{TT,n)/L'^~^~^ as a function of A. In short simulations, 
taking only a few days on an 800 MHz Pentium-Ill CPU, 
we find the quantum critical point at A = 0.396, which 
again compares well with earlier results [p^ . 

To summarize, we have presented Monte Carlo algo- 
rithms for the direct calculation of the free energy of 
a quantum system, based on a stochastic series expan- 
sion representation of the partition function. Our algo- 
rithms employ a variant of the Wang-Landau sampling 
to achieve a flat histogram and provide the free energy as 
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